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Abstract — We study the impact of sampling theorems on the 
fidelity of sparse image reconstruction on the sphere. We discuss 
how a reduction in the number of samples required to represent 
all information content of a band-limited signal acts to improve 
the fidelity of sparse image reconstruction, through both the 
dimensionality and sparsity of signals. To demonstrate this result 
we consider a simple inpainting problem on the sphere and 
consider images sparse in the magnitude of their gradient. We 
develop a framework for total variation (TV) inpainting on 
the sphere, including fast methods to render the inpainting 
problem computationally feasible at high-resolution. Recently a 
new sampling theorem on the sphere was developed, reducing the 
required number of samples by a factor of two for equiangular 
sampling schemes. Through numerical simulations we verify the 
enhanced fidelity of sparse image reconstruction due to the more 
efficient sampling of the sphere provided by the new sampling 
theorem. 

Index Terms — Spheres, harmonic analysis, sampling methods, 
compressive sensing. 

I. Introduction 

IMAGES are observed on a spherical manifold in many 
fields, from astrophysics (e.g. [ 1 ]) and biomedical imaging 
(e.g. (2)), to computer graphics (e.g. (3j) and beyond. In many 
of these settings inverse problems arise, where one seeks to 
recover an unknown image from linear measurements, which 
may be noisy, incomplete or acquired through a convolution 
process, for example. Such inverse problems are typically 
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solved by assuming some regularising prior on the unknown 
image to be recovered. 

Sparsity priors have received a lot of attention recently, 
since (a) they have been shown to be an effective and versatile 
approach for representing many real-world signals and (b) a 
sound theoretical foundation is provided by the emerging and 
rapidly evolving theory of compressive sensing (4J— (6). On the 
sphere, compressive sensing has been considered for signals 
sparse in the spherical harmonic domain [7 ], however a general 
theoretical framework does not yet exist for signals sparse in 
spatially localised representations. Nevertheless, sparse image 
reconstruction on the sphere in alternative representations, 
such as a set of overcomplete dictionaries, may still be 
considered; indeed, such an approach has been shown to be 
very effective (8). 

Although compressive sensing goes beyond Nyquist sam- 
pling, the Nyquist limit nevertheless defines the benchmark 
from which compressive sensing improvements are relative. 
Compressive sensing results are thus tightly coupled to the 
underlying sampling theorem on the manifold of interest. On 
the sphere, unlike Euclidean space, the number of samples 
required in the harmonic and spatial domains differ, with 
different sampling theorems on the sphere requiring a differ- 
ent number of samples in the spatial domain. Consequently, 
the sampling theorem adopted influences the performance of 
sparse signal reconstruction on the sphere. Studying the impact 
of sampling theorems on the sphere on the performance of 
sparse signal reconstruction is the focus of the current article. 

When considering signal priors that incorporate spatially 
localised information (for example directly in real space, in 
the magnitude of the gradient of signals, or through a wavelet 
basis or overcomplete dictionary), the sampling theorem that 
is adopted becomes increasingly important. Recently, a new 
sampling theorem on the sphere was developed by two of 
the authors of the current article for equiangular sampling 
schemes (9), reducing Nyquist sampling on the sphere by a 
factor of two compared to the canonical approach (10|, (TTJ. 
The reduction in the number of samples required to represent 
a band-limited signal on the sphere has important implications 
for sparse image reconstruction. 

To gain some intuition regarding these implications, we 
appeal to standard compressive sensing results in Euclidean 
space, where the ratio of the number of measurements M 
required to reconstruct a sparse image, to its dimensionality 
N, goes as M/N ex K |5j, (12), where K is the sparsity 
measure of the image (i.e. the number of non-zero coefficients 
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in some sparse representation)^] If one is not concerned with 
the number of measurements required to achieve a given re- 
construction fidelity but rather with the best fidelity for a given 
number of measurements, then this suggests reconstruction 
fidelity improves with decreasing dimensionality of the signal 
N and with decreasing sparsity K. 

Both of these quantities, dimensionality and sparsity, are 
related to the number of samples required to capture all 
information content of the underlying signal, as prescribed by 
the adopted sampling theorem. Spatial dimensionality is given 
identically by the number of samples of the sampling theorem. 
For any sparse representation of an image that captures spa- 
tially localised information, the sparsity of the signal is also 
directly related to spatial sampling. For example, in a wavelet 
representation, wavelets are located on each sample point. A 
less dense dictionary of wavelet atoms required to span the 
space will lead to a more sparse representation of images 
when the sparsity is computed in an analysis approach, i.e. 
as the number of non-zero projections of the signal onto the 
wavelet atoms. We concentrate on such analysis priors here, as 
suggested by the recent evolution of compressive sensing with 
redundant dictionaries fl2| . This argument can be extended to 
sparsity in the gradient and, in fact, all sparsity measures that 
capture spatially localised signal content. Consequently, for 
images sparse in a spatially localised representation, the ability 
to represent a band-limited signal on the sphere with fewer 
samples while still capturing all of its information content 
will improve the fidelity of sparse image reconstruction by 
enhancing both the dimensionality and sparsity of signals. 

In this article we study the implications of a new sampling 
theorem (9) for sparse image reconstruction on the sphere. 
We verify the hypothesis that a more efficient sampling of 
the sphere, as afforded by the new sampling theorem (9), 
enhances the fidelity of sparse image reconstruction through 
both the dimensionality and sparsity of signals. To demonstrate 
this result we consider a simple inpainting problem, where 
we recover an image on the sphere from incomplete spatial 
measurements. We consider images sparse in the magnitude 
of their gradient, as an illustration of the general setting, 
and develop a framework for total variation (TV) inpainting 
on the sphere. Solving these problems is computationally 
challenging; hence we develop fast methods for this purpose. 
Our framework is general and is trivially extended to other 
sparsity priors that incorporate spatially localised information. 

The remainder of the article is structured as follows. In 
Section [TT] we concisely review the harmonic structure of the 
sphere and corresponding sampling theorems. We develop a 
framework for TV inpainting on the sphere in Section [Till In 



showing the enhanced fidelity of sparse image reconstruction 
provided by a more efficient sampling of the sphere. Conclud- 



ing remarks are made in Section VI 



Section [IV] we describe algorithms for solving the optimisation 
problems on the sphere that arise in our TV inpainting 
framework. Numerical simulations are performed in Section [VJ 

typically the mutual coherence of the measurement and sparsifying 
operators also plays a role (5). However, in Euclidean space, as on the sphere, 
discrete inner products can be related to the (unique) continuous inner product 
(via a sampling theorem). Consequently, the measure of coherence is invariant 
to the choice of sampling theorem (the coherence is defined through the 
continuous inner product, which is the same for all sampling theorems). For 
the purpose of comparing sampling theorems on the sphere, we can thus safely 
neglect the impact of coherence. 



II. Sampling on the Sphere 

A sampling theorem on the sphere states that all information 
in a (continuous) band-limited signal is captured in a finite 
number of samples in the spatial domain. Since a (continuous) 
band-limited signal on the sphere may be represented by a 
finite harmonic expansion, a sampling theorem on the sphere 
is equivalent to an exact prescription for computing a spherical 
harmonic transform from a finite set of spatial samples. In this 
section we review the harmonic structure of the sphere, before 
discussing sampling theorems on the sphere. 

A. Harmonic structure of the sphere 

We consider the space of square integrable functions on the 
sphere L 2 (S 2 ), with the inner product of x,y G L 2 (S 2 ) defined 

(x, y)= / dSl(O,<p)x(0,<p)y*(O,<p), 
Js 2 

where dQ(9,ip) = sinOdOdip is the usual invariant measure 
on the sphere and (9,ip) denote spherical coordinates with 
colatitude 9 G [0, ir] and longitude ip G [0, 27r). Complex 
conjugation is denoted by the superscript *. The canonical 
basis for the space of square integrable functions on the sphere 
is given by the spherical harmonics Y^ m G L 2 (S 2 ), with 
natural i G N, integer m G Z and \m\ < £. Due to the 
orthogonality and completeness of the spherical harmonics, 
any square integrable function on the sphere x G L 2 (S 2 ) may 
be represented by its spherical harmonic expansion 



(1) 



£=0 m=-£ 



where the spherical harmonic coefficients are given by the 
usual projection onto each basis function: 

X£m = (X, Yim) = / dQ(9, <p) x(0, if) Yg m {9 , if) . 

Js 2 

Throughout, we consider signals on the sphere band-limited 
at L, that is signals such that X£ m = 0, W > L, in which 
case the summation over i in ([I]) may be truncated to the 
first L terms. Finally, note that the harmonic coefficients of 
a real function on the sphere satisfy the conjugate symmetry 
relation x* im = (-l) m X£ —m, which follows directly from the 
conjugate symmetry of the spherical harmonics. 



B. Sampling theorems on the sphere 

Sampling theorems on the sphere describe how to sample a 
band-limited signal x so that all information is contained in a 
finite number of samples N . Moreover, a sampling theorem on 
the sphere effectively encodes an exact quadrature rule for the 
integration of band-limited functions | 9 ], [10]. We denote the 
concatenated vector of N spatial measurements by x e C N 
and the concatenated vector of L 2 harmonic coefficients by 

T 2 

x G C . The number of spatial and harmonic elements, N 
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and L 2 respectively, may differ (and in fact do differ for all 
known sampling theorems on the sphere). 

Before discussing different sampling theorems on the 
sphere, we define a generic notation to describe the harmonic 
transform corresponding to a given sampling theorem. A 
sampling theorem describes how to compute the spherical 
harmonic transform of a signal exactly. Since the spheri- 
cal harmonic transform and inverse are linear, we represent 
the forward and inverse transform by the matrix operators 
reC L 2 xAr and A e £NxL 2 respectively The spherical 

harmonic coefficients of a sampled signal (i.e. image) on the 
sphere x are given by the forward transform 

x = Yx . 

while the original signal is recovered from its harmonic 
coefficients by the inverse transform 

x = Ax . 

Different sampling theorems then differ in the definition of 
A, T and the number of spatial samples N. By definition, all 
sampling theorems give exact spherical harmonic transforms, 
implying TA = t L 2, where 1^ is the k x k identity matrix. 
However, for all sampling theorems on the sphere the number 
of samples required in the spatial domain exceeds the number 
of coefficients in the harmonic domain (i.e. N > L 2 ), hence 
Ar 7^ tjsf. Consequently, for the N sample positions of a 
sampling theorem, an arbitrary set of sample values does 
not necessarily define a band-limited signal (contrast this to 
the discrete Euclidean setting where a finite set of samples 
uniquely defines a band-limited signal). Note also that the 
adjoint inverse (forward) spherical harmonic transform differs 
to the forward (inverse) spherical harmonic transform in the 
discrete setting. 

For an equiangular sampling of the sphere, the Driscoll & 
Healy (DH) (10| sampling theorem has become the standard, 
requiring TVdh = 2L(2L — 1) ~ 4L 2 samples on the sphere 
to represent exactly a signal band-limited in its spherical 
harmonic decomposition at L. Recently, a new sampling 
theorem for equiangular sampling schemes has been de- 
veloped by McEwen & Wiaux (MW) (9), requiring only 
]Vmw = (L — 1)(2L — 1) + 1 ~ 2L 2 samples to represent a 
band-limited signal exactly. No sampling theorem on the 
sphere reaches the optimal number of samples suggested by 
the L 2 dimension of a band-limited signal in harmonic space 
(although the MW sampling theorem comes closest to this 
bound). The MW sampling theorem therefore achieves a more 
efficient sampling of the sphere, with a reduction by a factor 
of approximately two in the number of samples required to 
represent a band-limited signal on the sphere^] 

2 Gauss-Legendre (GL) quadrature can also be used to construct an efficient 
sampling theorem on the sphere, with Nql — L(2L — 1) ~ 2L 2 samples 
(see e.g. (9]). The MW sampling theorem nevertheless remains more efficient, 
especially at low band-limits. Furthermore, it is not so straightforward to 
define the TV norm on the GL grid since it is not equiangular. Finally, algo- 
rithms implementing the GL sampling theorem have been shown to be limited 
to lower band-limits and less accurate than the algorithms implementing the 
MW sampling theorem |9 |. Thus, we focus on equiangular sampling theorems 
only in this article. 



Fast algorithms have been developed to compute forward 
and inverse spherical harmonic transforms rapidly for both 
the DH (10), (TTJ and MW (9) sampling theorems. These 
fast algorithms are implemented, respectively, in the publicly 
available SpharmonicKit£ package and the Spin Spherical 
Harmonic Transform (SSHT^J package and are essential to 
facilitate the application of these sampling theorems at high 
band-limits. 



III. Sparse Image Reconstruction on the Sphere 

A more efficient sampling of a band-limited signal on the 
sphere, as afforded by the MW sampling theorem, improves 
the quality of sparse image reconstruction for images that are 
sparse in a spatially localised measure. To demonstrate this 
result we consider a simple inpainting problem on the sphere 
and consider images sparse in the magnitude of their gradient. 
We develop a framework for total variation (TV) inpainting 
on the sphere, which relies on a sampling theorem and its 
associated quadrature rule to define a discrete TV norm on the 
sphere. Firstly, we define the discrete TV norm on the sphere, 
before secondly defining finite difference gradient operators 
on the sphere. Thirdly, we discuss the TV inpainting problem. 



A. TV norm on the sphere 

We define the discrete TV norm on the sphere by 



N e -1 N^-l 

\\x\\tv = £ E l(0 t ) \Vx\ , 

t=0 p=0 



(2) 



where t and p index the equiangular samples in 6 and (p 
respectively, with the number of samples associated with a 
given sampling theorem denoted in each dimension by Nq 
and respectively. The discrete magnitude of the gradient 
is defined by 



(S e x) 



sin z 0, 



(S^x)' 



(3) 



to approximate the continuous magnitude of the gradient 




by finite differences. The finite difference operators 5q and 5<p 
are defined explicitly in the following subsection. The con- 
tribution to the TV norm from the magnitude of the gradient 
for each pixel value is weighted by the quadrature weights 
q(0t) of the sampling theorem adopted in order to approximate 



http ://www.cs. dartmouth . edu/~ geelong/sphere/ 
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continuous integration]^] The inclusion of the weights q(0 t ) 
also regularises the sin 6 term that arises from the definition of 
the gradient on the sphere, eliminating numerical instabilities 
that this would otherwise cause. 

B. Gradient operators on the sphere 

The finite difference operators Sq and 5^ defined on the 
sphere appear in the definition of the discrete magnitude of 
the gradient given by and thus are required to compute the 
discrete TV norm on the sphere. Furthermore, as we shall see, 
to solve the TV inpainting problems outlined in the following 
subsection, the adjoints of these operators are also required. 
We define these operators and adjoints explicitly here. 

The operator 5q is defined sample-wise by 




■x tiP , * = 0,1,--- ,N e -2 and Vp 
t = N e - 1 and \/p 



with adjoint 



t, P 



u t -i, p - 
u t-l,p, 



t = and Vp 
Ut lP , t = 1, • • • Nq — 2 and \/p . 
t = N e - 1 and Mp 



Note that this definition is identical to the typical definition of 
the corresponding operator on the plane |T3| . The operator 5^ 
is defined sample-wise by 




with adjoint 

(*>kp = 



p 
p 



Vt,p-1 - v tiP , 



0,1,- •• -2 and Vt 
N^ — 1 and Vt 



p 
p 



and \ft 



1 and \/t 



Since the sphere is periodic in ip, we define the corresponding 
finite difference operator to also be periodic. The finite differ- 
ence operator and adjoint in ip therefore differ to the typical 
definition on the plane fT3| . 

The TV norm on the sphere may then be seen as the sum 
of the magnitude of the weighted gradient 



where 



\ X \\TV 



Nq — 1 n^-i 
p=0 



£ 

t=o 



t, P \ 



t, P 



= ( U t,p ' 



t,p 



1/2 



5 If the band-limiting operator T = Ar £ C NxN were applied to |Vsc| 
in j2j, then the finite summation of (|2j would give an exact quadrature for 
the integral of the continuous function underlying the associated samples of 
the band-limited |Vsc|. However, introducing the operator T in J2J would 
make solving the optimisation problems defined subsequently problematic 
and would also prohibit passing the quadrature weights inside the gradient 
to eliminate numerical instabilities due to the sin# term. Consequently, we 
adopt the definition of the discrete TV norm on the sphere given by |2j. 
In any case, numerical experiments have shown that ||sc||tv i s identical 
to E^ P Vas | to m achine precision for the particular test images 

considered in Section |V- A | Thus, the discrete TV norm defined by |2j can be 
though of as an accurate proxy for J g2 dfl Y|Vsc|. 



for 



Vx . 

The weighted gradient operator is defined by 



V 



where the weighted finite difference operators are defined by 



{6e) t =q{0 t )(5 e ) 



and 



q(Qt) 

sin t 



(5. 



Notice how the inclusion of the weights q(0 t ) regularises the 
sin^ term that arises from the definition of the gradient on 
the sphere, eliminating numerical instabilities that this would 
otherwise cause. If 6 t = 7r, corresponding to the South pole of 
the sphere, then (S ip x) tp = and thus we define (5<px) t = 
to avoid dividing by sin 6^ = 0. Note that the MW sampling 
theorem includes a sample on the South pole, while the DH 
sampling theorem does not (neither sampling theorem includes 
a sample on the North pole). The adjoint weighted gradient 
operator is then applied as 



where the adjoint operators 5 \ and 5^ follow trivially from S J 



and 5^. 



C. TV inpainting on the sphere 

We consider the measurement equation 

y — + n , 

where M noisy real measurements y G R M of the underlying 
real image on the sphere x G R N are made. The matrix 
implementing the measurement operator G R MxN repre- 
sents a uniformly random masking of the image, with one 
non-zero, unit value on each row specifying the location of 
the measured datum. The noise n G R M is assumed to be 
independent and identically distributed (iid) Gaussian noise, 
with zero mean and variance o\. We assume that the image 
x is sparse in the norm of its gradient and thus attempt to 
recover x from measurements y by solving the following TV 
inpainting problem directly on the sphere: 

x* = argmin ||#||tv such that \\y — <&x\\2 < e . (4) 
x 

The square of the residual noise follows a scaled x 2 distribu- 
tion with M degrees of freedom, i.e. \\y—<frx*\\2 ~ cr^x 2 (M). 
Consequently, we choose e 2 to correspond to the (100a)th 
percentile of this distribution, giving a probability a that pure 
noise produces a residual noise equal to or smaller than the 
observed residual. Note that the data constraint in ^ is given 
by the usual £2 -norm, which is appropriate for Gaussian noise 
on a discrete set of measurements. Although we consider band- 
limited signals, we have not imposed this constraint when 
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solving Q. Consequently, x* will not necessarily be band- 
limited at L and we impose this constraint on the solution 
by performing a forward and inverse spherical harmonic 
transform: x\ = Tx*, where the band-limiting operator is 
defined by T = Ar G C NxN . 

As discussed already, for images sparse in a measure that 
captures spatially localised information, such as the TV norm, 
a more efficient sampling of the signal enhances sparsity. 
Furthermore, when recovering signals in the spatial domain 
directly, the dimensionality of the signal is also enhanced 
by a more efficient sampling. These two effects both act to 
improve the fidelity of sparse image reconstruction. Thus, 
the more efficient sampling of the MW sampling theorem 
when compared to the DH sampling theorem will improve 
the fidelity of sparse image reconstruction when solving the 
TV inpainting problem given by ([4]). We verify these claims 
with numerical experiments in Section [V] 

No sampling theorem on the sphere reaches the optimal 
number of samples in the spatial domain suggested by the 
L 2 dimensionality of the signal in the harmonic domain. We 
may therefore optimise the dimensionality of the signal that 
we attempt to recover by recovering its harmonic coefficients 
x directly. We do so by solving the following TV inpainting 
problem in harmonic space: 

x'* = argmin ||A'£'|| T v such that \\y - $A'x\\ 2 < e . (5) 
x 

We impose reality of the recovered signal by explicitly im- 
posing conjugate symmetry in harmonic space through the 
conjugate symmetry extension operator II G C L xL ( L + 1 )/ 2 9 
where A' = All. The full set of harmonic coefficients of x are 
given by x = Ux , where x G C L( ^ L+1 ^ 2 are the harmonic 
coefficients for the spherical harmonic azimuthal index m non- 
negative only. The image on the sphere is then recovered from 
its harmonic coefficients by x* = A'x'*. By solving the TV 
inpainting problem directly in harmonic space, we naturally 
recover a signal band-limited at L. 

When solving the TV inpainting problem ([5]) directly in 
harmonic space, the dimensionality of the recovered signal is 
optimal and identical for both sampling theorems. However, 
the sparsity of the signal with respect to the TV norm remains 
enhanced for the MW sampling theorem when compared to 
the DH sampling theorem. Consequently, the MW sampling 
theorem will improve the fidelity of sparse image reconstruc- 
tion when solving the TV inpainting problem given by ([5]), 
although through sparsity only and not also dimensionality. We 
verify these claims with numerical experiments in Section |V| 

Note that if a band-limit constraint were explicitly imposed 
in problem ([?]), then the two problems would be equivalent, 
however, this would involve applying the band-limiting op- 
erator T = AT, complicating the problem and increasing 
the computational cost of finding a solution, while providing 
no improvement over (J5J. In the current formulation of these 
two optimisation problems, problem ^ has the advantage of 
simplicity, while problem §5§ is the simplest formulation that 
optimises dimensionality. 



IV. Algorithms 

We solve the TV inpainting problems on the sphere given by 
^ and §5§ using iterative convex optimisation methods. Solv- 
ing the TV inpainting problem in harmonic space poses two 
challenges as we go to high band-limits (i.e. high-resolution). 
Firstly, we require as an input to the convex optimisation 
algorithm an upper bound on the inverse transform operator 
norm, which is challenging to compute at high-resolution. 
We describe a method to compute the operator norm at 
high-resolution, which, crucially, does not require an explicit 
computation of A. Secondly, the inverse spherical harmonic 
transform A and its adjoint operator A^ must be applied 
repeatedly in the iterative algorithm. Fast algorithms are 
essential to perform forward and inverse spherical harmonic 
transforms at high-resolution and have been developed for 
both the DH (TO), (TTJ and MW |9] sampling theorems. 
To solve the inpainting problem at high-resolution we also 
require a fast adjoint inverse transform. We thus develop 
fast algorithms to perform the adjoint forward and adjoint 
inverse spherical harmonic transforms corresponding to the 
MW sampling theorem. Since we predict the MW sampling 
theorem to be superior to the DH sampling theorem for 
sparse image reconstruction on the sphere (a prediction that is 
validated by numerical experiments performed at low band- 
limits (i.e. low-resolution) in Section [V]), we develop fast 
adjoint algorithms for the MW sampling theorem only. These 
methods then render the computation of solutions to the TV 
inpainting problems feasible at high-resolution for the MW 
sampling theorem. 



A. Convex optimisation 

We apply the Douglas-Rachford proximal splitting algo- 
rithm [14] to solve the convex optimisation problems Q and 
(5)|^We describe only how to solve problem §5§ as problem 
( 4 ) can be solved in the same way (by replacing A' with the 
identity matrix In and by replacing x with x). 

The Douglas-Rachford algorithm fl4| is based on a splitting 
approach that requires the computation of two proximity 
operators (HJ. In our case, one proximity operator is based 
on the TV norm ||A' • ||tv an d the other on the data constraint 
||y-$A'.|| 2 <€. 

In the case of an image on the plane, the proximity operator 
based on the TV norm may be computed using, for example, 
the method described in |T3| or in |T6| . For an image on the 
sphere, the same methods can be used after introducing the 
following modifications. In |T6| the algorithm to compute the 
proximity operator of the TV norm is described in terms of a 
linear operator £, its adjoint C\ and two projections onto a 
set V and a set C. In our case, the linear operator C and its 
adjoint & may be redefined as 



-A'+vt 



6 We use Douglas-Rachford splitting since this does not require differentia- 
bility of the objective function and allows us to solve constrained optimisation 
problems where we adopt an indicator function to represent the measurement 
constraint. 
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and 

where the set V is the set of weighted gradient-pairs (u,v) 
such that u 2 p + < 1 and C is simply given by the space 
of the recovered vector x. 

The second proximity operator, related to the data constraint 
1 1 y — <$>A f '|| 2 < e, is computed using the method described 
in p7) directly. 

5. Operator norm bound 

The convex optimisation algorithm requires as input upper 
bounds for the norms of the operators that appear in the 
problem. The calculation of these norms is in most cases 
straightforward, however the calculation of the inverse spher- 
ical harmonic transform operator norm, defined by 

||A|| 2 = max ||A4|| 2 , 
l|£|| 2 =i 

can prove problematic. At low-resolution ||A||2 may be com- 
puted explicitly, however this is not feasible at high-resolution 
since even computing and storing A explicitly is challenging. 

We develop a method here to estimate this norm for the MW 
sampling theorem without computing A explicitly. We seek 
a sampled function on the sphere x = Ax that maximises 
1 1 x 1 1 2, while satisfying the constraint ||£||2 = 1. By the 
Parseval relation and the sampling theorem on the sphere, this 
constraint may be rewritten: 

||£|| 2 = 1 => (x,x) = l 

Parseval 

Sampling theorem 

where x u G R Nu contains samples of x, sampled at a 
resolution sufficient to represent x 2 , i.e. corresponding to 
band-limit 2L — 1 (so that an exact quadrature may be used to 
evaluate (x, x) from a discrete set of samples), Q u G R N " xN * 
is the matrix with corresponding quadrature weights along its 
diagonal, and where N u ~ 2(2L — l) 2 . Since we know that the 
quadrature weights for the MW sampling theorem are closely 
approximated by sin# (9J, the signal that maximises ||aj||2 
while satisfying the constraint x u ^Q u x u = 1 has its energy 
centred as much as possible on the South pole since this is 
where the quadrature weights are smallest (recall that the MW 
sampling scheme does not contain a sample on the North pole). 
This signal is given by the band-limited Dirac delta function 
centred on the South pole (see e.g. fT8| for the definition of the 
band-limited Dirac delta function on the sphere). The spherical 
harmonic coefficients of this band-limited Dirac delta function 
5 L G L 2 (S 2 ) are given by 

where k is a normalisation factor chosen to ensure ||(S L ||2 = 1 
and Sij is the Kronecker delta symbol. The norm of the inverse 
spherical harmonic transform operator may then be computed 



80 




L 



Fig. 1. Explicit calculation of the inverse spherical harmonic transform 
operator norm ||A||2 and estimation by the method outlined in the text, at 
low-resolution. The solid red line shows the estimated norm for all band- 
limits L, while the solid blue circles show the values computed explicitly for 
L £ {2, 4, 8, 16, 32}. The estimated norm agrees with the actual norm very 
well. 

by || A || 2 — ||A^ L ||2, which, crucially, does not require an 
explicit computation of A, merely its application. 

In Figure [T] we compute ||A||2 by the method outlined here 
and from A explicitly, for low-resolution. We find that the 
method to estimate the norm of the inverse spherical harmonic 
transform operator outlined here estimates the actual norm 
very well. 

We also derived an upper bound for the norm of this 
operator for the MW sampling theorem. However, the bound 
we derived is not tight and we found empirically that the 
method outlined here to estimate the norm itself, rather than 
a bound, is very accurate and improved the performance of 
the optimisation algorithm considerably when compared to 
a non- tight bound. Although we do not prove so explicitly, 
we conjecture that the method outlined here gives the inverse 
transform operator norm exactly. 

C. Fast adjoint spherical harmonic transforms 

Standard convex optimisation methods require not only the 
application of the operators that appear in the optimisation 
problem but often also their adjoints. Moreover, these methods 
are typically iterative, necessitating repeated application of 
each operator and its adjoint. Thus, to solve optimisation 
problems that incorporate harmonic transform operators, like 
the harmonic space TV inpainting problem given by ([5]), fast 
algorithms to apply both the operator and its adjoint are 
required to render high-resolution problems computationally 
feasible. 

Here we develop fast algorithms to perform adjoint forward 
and adjoint inverse spherical harmonic transforms for the 
MW sampling theorem. Although we only require the adjoint 
inverse transform in this article, for the sake of completeness 
we also derive a fast adjoint forward transform. Similarly, 
although we only consider scalar functions in this article, for 
the sake of completeness we derive fast adjoint algorithms 
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for the spin setting. A spin function on the sphere transforms 



as s f r (0,Lp) 



3 f(0,(f) under a local rotation by 



X £ [0, 27r), where the prime denotes the rotated function. 
It is important to note that the rotation considered here is not 
a global rotation on the sphere but rather a rotation by \ m me 
tangent plane at (0, (p) (see e.g. |9J for further details). In the 
expressions for the fast algorithms derived below, the standard 
scalar case follows simply by setting s = 0. These fast adjoint 
algorithms are implemented in the publicly available SSHtQ 
package (9). 

The fast adjoint inverse spherical harmonic transform for 
the MW sampling theorem follows by taking the adjoint of 
each stage of the fast inverse transform [ 9 ] and applying these 
in reverse order. The final stage of the fast inverse transform 
involves discarding out-of-domain samples and has adjoint 



a f(0t,<Pp), *G{0,1,...,L-1} 



0. 



t e {L, ...,2L-2} 

The second stage of the fast adjoint inverse transform is given 
by 

2L-2 2L-2 

t=o p=0 

which may be computed rapidly using fast Fourier transforms 
(FFTs). The final stage of the fast adjoint inverse transform is 
given by 



47T 



L-1 



/ j ^m'm <-*m' ,-s s 1 mm' •> 
m'=-(L-l) 

where A^ n = d^ n (7r/2) are the Wigner d- functions evalu- 
ated for argument 7r/2 (see e.g. |19|). This final calculation 
dominates the overall asymptotic complexity of the fast adjoint 
inverse transform, resulting in an algorithm with complexity 
G(L 3 ). 

The fast adjoint forward spherical harmonic transform for 
the MW sampling theorem follows by taking the adjoint of 
each stage of the fast forward transform [ 9 ] and applying these 
in reverse order. The first stage of the fast adjoint forward 
transform is given by 

s Ow f =(-l)*r(™+*) 

£=0 

The next stage is given by the (reflected) convolution 

L-1 

S ^mm" t =27T ^ S G W f W \w! ~ TU' ' ') , 

m'=-(L-l) 

which is self-adjoint, followed by the inverse Fourier transform 
in 

1 L_1 

F USA — - F z 1 " e irn ' 9t 

s ± m \vt) — 2^ 2 / _j s 1 mm' c i 

m / = -(L-l) 



which may be computed rapidly using FFTs. The next stage 
consists of the adjoint of the periodic extension of a function 
on the sphere performed in the forward transform and is given 
by 

8 Fj(e t ) = 

+ (-l) m+s a Fj{6 2L _ 2 _ t ) , t e {0, 1, . . . , L - 2} . 
,Fj(0 t ), t = L-l 

The final stage consists of the Fourier transform in ip 

L-1 



irriLpp 



http://www.jasonmcewen.org/ 



m=-(L-l) 

which may be computed rapidly using FFTs. The first cal- 
culation dominates the overall asymptotic complexity of the 
fast adjoint forward transform, resulting in an algorithm with 
complexity 0(L 3 ). 

V. Simulations 

We perform numerical experiments to examine the impact 
of a more efficient sampling of the sphere when solving the 
TV inpainting problems defined in Section [TTTJ Firstly, we 
perform a low-resolution comparison of reconstruction fidelity 
when adopting the DH and MW sampling theorems, where the 
predicted improvements in reconstruction fidelity provided by 
the MW sampling theorem are verified in practice. Secondly, 
we perform a single simulation to illustrate TV inpainting at 
high-resolution on a realistic test image. 

A. Low -re solution comparison on band-limited images 

A test image is constructed from Earth topography data. 
The original Earth topography data are taken from the Earth 
Gravitational Model (EGM2008) publicly released by the 
U.S. National Geospatial-Intelligence Agency (NGA) EGM 
Development Team|j To create a band-limited test signal 
sparse in its gradient, the original data are thresholded at their 
midpoint to create a binary Earth map (scaled to contain zero 
and unit values), which is then smoothed by multiplication in 
harmonic space with the Gaussian G^ m = exp(— £ 2 a s ), with 
cr s = 0.002, to give a signal band-limited at L = 32. The 
resulting test image is displayed in Figure [2] Let us stress that 
this test image is constructed to satisfy the assumptions of our 
theoretical framework, i.e. the case of band-limited images that 
are sparse in their gradient. This is necessary to evaluate the 
theoretical predictions based on our framework. A realistic test 
image is considered in the following subsection. 

Measurements of the test image are taken at uniformly 
random locations on the sphere, as described by the mea- 
surement operator <£, in the presence of Gaussian iid noise 
with standard deviation a n = 0.01. Reconstructed images on 
the sphere are recovered by solving the inpainting problems 
in the spatial and harmonic domains, through ^ and §5§ 
respectively, using both the DH and MW sampling theorems, 

8 These data were downloaded and extracted using the tools available from 
Frederik Simons' webpage: http://www.frederik.net 



8 



IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 22, NO. 6, JUNE 2013 




Fig. 2. Test image of Earth topographic data constructed to be sparse in its 10 "/ 

gradient and band-limited at L = 32. This image constitutes the ground truth > y 
in our numerical experiments. Here and subsequently data on the sphere are 5 

displayed using the Mollweide projection, with zero values shown in black, ; ■ '■ ■ ■ j 

unit values shown in yellow, and the colour of intermediate values interpolated 0.25 0.5 0.75 1 1 .25 1 .5 1 .75 2 

between these extremes. M / L 



giving four reconstruction techniques. The bound e is deter- 
mined from a — 0.99. We consider the measurement ratios 
M/L 2 e {1/4, 1/2, 1, 3/2, Nuw/L 2 ~ 2} (recall that L 2 
is the dimensionality of the signal in harmonic space). The 
measurement ratio M/L 2 = Nmw/L 2 ~ 2 corresponds to 
complete coverage for the MW sampling theorem, i.e. Nyquist 
rate sampling on the MW grid. 

Typical reconstructed images are shown in Figure [3] for 
the four reconstruction techniques. For each reconstruction 
technique and measurement ratio M/L 2 , we perform ten 
simulations for random measurement operators and noise. To 
quantify the error of reconstruction, we compute the signal-to- 
noise-ratio SNR = 201og(||£|| 2 /||£* - x\\ 2 ) (defined in har- 
monic space to avoid differences due to the number of samples 
of each sampling theorem). Note that the standard ^ 2 -norm is 
used in the definition of the SNR given the discrete nature of 
harmonic space on the sphere. Reconstruction performance, 
averaged over these ten simulations, is shown in Figure [4] 

When solving the inpainting problem in the spatial domain 
through ^ we see a large improvement in reconstruction 
quality for the MW sampling theorem when compared to 
the DH sampling theorem. This is due to the enhancement 
in both dimensionality and sparsity afforded by the MW 
sampling theorem in this setting. When solving the inpaint- 
ing problem in the harmonic domain through ([5]) we see a 
considerable improvement in reconstruction quality for each 
sampling theorem, since we optimise the dimensionality of 
the recovered signal by going to harmonic space. For harmonic 
reconstructions, the MW sampling theorem remains superior to 
the DH sampling theorem due to the enhancement in sparsity 
(but not dimensionality) that it affords in this setting. All of 



the predictions made in Section [Hi] are thus exhibited in the 
numerical experiments performed in this section. In all cases, 
the superior performance of the MW sampling theorem is 
clear. 



B. High-resolution illustration on a realistic image 

In this section we perform a single simulation to illustrate 
TV inpainting at high resolution. Furthermore, we also con- 
sider a more realistic test image. Since we develop fast adjoint 
algorithms for the MW sampling theorem only (due to its 



Fig. 4. Reconstruction performance for the DH (green/diamonds) and MW 
(red/circles) sampling theorems, when solving the TV inpainting problem in 
the spatial (dot-dashed line) and harmonic domain (solid line). The MW 
sampling theorem provides enhancements in reconstruction quality when 
compared to the DH sampling theorem, due to dimensionality and sparsity 
improvements in spatial reconstructions, and due to sparsity (but not dimen- 
sionality) improvements in harmonic reconstructions. 



superiority), we therefore use only the MW sampling theorem 
for the high-resolution inpainting simulation performed here. 
A high-resolution test image is constructed from the same 



Earth topography data described in Section |V-A| Since the 
original data are defined in harmonic space, we first simply 
truncate the harmonic coefficients to yield an image band- 
limited at L = 128. In practice, acquired images may not 
necessarily be band-limited. We thus process the data to 
construct a test image that is not band-limited. We construct 
such a test image defined by height above sea-level (i.e. we 
threshold all oceans and trenches, while the continents and 
mountains remain unaltered). The abrupt transition between 
the oceans and the continents results in an image that is indeed 
not band- limited (as verified numerically). Furthermore, the 
continents and mountainous regions result in a test image that 
is not highly sparse in its gradient. The resulting realistic test 
image is shown in Figure [5] (a). 

The same measurement procedure as outlined previously 
is applied to take noisy, incomplete measurements of the 
data for a range of measurement ratios M/L 2 . The inpainted 
images are recovered by solving the inpainting problem in 
harmonic space through §5§ using the MW sampling theorem. 
To solve the inpainting problem for these high-resolution 
simulations we use the estimator of the inverse transform 



norm ||A|| 2 described in Section IV-B and the fast adjoint 



harmonic transform algorithms defined in Section IV-C Using 
these fast algorithms, combined with recent optimisations of 
the SSHT package, it takes approximately 10 minutes to solve 
the inpainting problem in harmonic space at L = 128 on a 
standard laptop (with a 1.8 GHz Intel Core i7 processor and 
4 GB of RAM). 

The inpainted images are shown in Figure [5] Since the 
original realistic test image is not band-limited, the previous 
SNR measure (which is defined in harmonic space to avoid 
a dependence on the number of samples of each sampling 
theorem) is not a meaningful error metric (computation of 
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(a) DH spatial for M/L 2 = 1/4 (b) DH harmonic for M/L 2 = 1/4 (c) MW spatial for M/L 2 = 1/4 (d) MW harmonic for M/L 2 = 1/4 




(e) DH spatial for M/L 2 = 1/2 (f) DH harmonic for M/L 2 = 1/2 (g) MW spatial for M/L 2 = 1/2 (h) MW harmonic for M/L 2 = 1/2 




(i) DH spatial for M/L 2 = 1 (j) DH harmonic for M/L 2 = 1 (k) MW spatial for M/L 2 = 1 (1) MW harmonic for M/L 2 = 1 




(m) DH spatial for M/L 2 = 3/2 (n) DH harmonic for M/L 2 = 3/2 (o) MW spatial for M/L 2 = 3/2 (p) MW harmonic for M/L 2 =3/2 




(q) DH spatial for M/L 2 ~ 2 (r) DH harmonic for M/L 2 ~ 2 (s) MW spatial for M/L 2 ~ 2 (t) MW harmonic for M/L 2 ~ 2 

Fig. 3. Inpainted images on the sphere recovered by solving the TV inpainting problems for a range of measurement ratios M/L 2 . The ground truth image 
is shown in Figure [2] The first and second columns of panels show the inpainted images recovered using the DH sampling theorem, while the third and fourth 
columns show the inpainted images recovered using the MW sampling theorem. The first and third columns of panels show inpainted images recovered by 
solving the inpainting problem in the spatial domain, while the second and fourth columns show images recovered by solving the inpainting problem in the 
harmonic domain. The final row of panels corresponds to measurement ratio M/L 2 = Nmw /L 2 ~ 2. The quality enhancements due to the MW sampling 
theorem and by solving the inpainting problem in harmonic space are both clear. 



the harmonic transform would indeed be affected by un- 
controlled aliasing). Instead, we use the analogous SNR 
measure defined in image space on the sphere, given by 
SNRi = 101og(ajtQaj/((^ - x)^Q(x* - x))), where we re- 
call Q is the matrix with quadrature weights on its diagonal. 
Just as for the definition of the TV norm, the inclusion of the 
weights for signals that are not band-limited provides only an 
intuitive approximation to continuous integration. The SNRi 
values for each inpainted image are displayed in Figure [5] 
from which it is apparent that SNRi increases with increasing 
measurement ratio. Moreover, it is clearly apparent by eye 
that our TV inpainting framework is effective when applied 
to realistic images that are not highly sparse in their gradient 
and that are not band-limited. 



VI. Conclusions 

The MW sampling theorem, developed only recently, 
achieves a more efficient sampling of the sphere than the 
standard DH sampling theorem: without any loss to the 
information content of the sampled signal, the MW sampling 
theorem reduces the number of samples required to represent 
a band-limited signal by a factor of two for an equiangular 
sampling. For signals sparse in a spatially localised measure, 
such as in a wavelet basis, overcomplete dictionary, or in 
the magnitude of their gradient, for example, a more efficient 
sampling enhances the fidelity of sparse image reconstruction 
through both dimensionality and sparsity. When a signal is 
recovered directly in the spatial domain, the MW sampling 
theorem provides enhancements in both dimensionality and 
sparsity when compared to the DH sampling theorem. By 
recovering the signal directly in harmonic space it is possible 
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(a) Ground truth 





(b) Inpainted image for M/L 2 = 1/4 (SNRi = 20.0dB) 
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(c) Inpainted image for M/L 2 = 1/2 (SNRi = 27.8dB) (d) Inpainted image for M/L 2 = 1 (SNRi = 37.0dB) 
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(e) Inpainted image for M/L 2 = 3/2 (SNRi = 38.5dB) (f) Inpainted image for M/L 2 - 2 (SNRi = 53.2dB) 

Fig. 5. Inpainting illustration for a realistic image at high-resolution (L = 128). The inpainted images are recovered by solving the inpainting problem in 
harmonic space using the MW sampling theorem for a range of measurement ratios M/L 2 . The SNRi of each recovered image is also displayed. 



to optimise its dimensionality, in which case the MW sampling 
theorem still provides an enhancement in sparsity but not in 
dimensionality. 

We verified these statements through a simple inpainting 
problem on the sphere, where we considered images sparse in 
their gradient. We built a framework and fast methods for total 
variation (TV) inpainting on the sphere. Using this framework 
we performed numerical experiments which confirmed our 
predictions: in all cases, the more efficient sampling provided 
by the MW sampling theorem improved the fidelity of sparse 
image reconstruction on the sphere. 
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